#!/usr/bin/env julia
"""
Simple diagnostic file generation test
"""

using Printf
using Dates
using Random

# Load modules directly
include(joinpath(@__DIR__, "..", "src", "Diagnostics", "FortranBinaryIO.jl"))
include(joinpath(@__DIR__, "..", "src", "Diagnostics", "DiagnosticFileFormat.jl"))
include(joinpath(@__DIR__, "..", "src", "Diagnostics", "DiagnosticWriters.jl"))

import .FortranBinaryIO as FBI
import .DiagnosticFileFormat as DFF
import .DiagnosticWriters as DW

println("="^80)
println("Testing GSI Diagnostic File Generation")
println("="^80)

# Create output directory
output_dir = joinpath(@__DIR__, "..", "results", "diagnostics")
mkpath(output_dir)

# Create some test observations
Random.seed!(42)
observations = DFF.ObservationRecord[]

for i in 1:100
    obs = DFF.ObservationRecord()
    obs.lat = rand() * 180.0 - 90.0
    obs.lon = rand() * 360.0
    obs.pressure = rand() * 1000.0 + 10.0
    obs.obs_value = 1013.0 + randn() * 5.0
    obs.obs_error = 1.0
    obs.background_value = obs.obs_value + randn() * 2.0
    obs.analysis_value = obs.obs_value + randn() * 1.0
    obs.ob_minus_background = obs.obs_value - obs.background_value
    obs.ob_minus_analysis = obs.obs_value - obs.analysis_value
    obs.qc_flag = Int32(0)
    obs.use_flag = Int32(1)
    obs.obs_type = Int32(120)
    obs.obs_subtype = Int32(181)
    obs.station_id = @sprintf("STN%05d", i)
    obs.time_offset = 0.0
    obs.inverse_obs_error = 1.0 / obs.obs_error
    obs.variational_qc_weight = 1.0

    push!(observations, obs)
end

println("Created $(length(observations)) observations")

# Create diagnostic output container
analysis_time = DateTime(2018, 8, 12, 12, 0, 0)
grid_size = (95, 57, 32)

diag = DW.DiagnosticOutput(output_dir, analysis_time, grid_size)
diag.observations["ps"] = observations

# Compute statistics
stats = DFF.compute_innovation_statistics(observations)
diag.innovation_stats["ps"] = stats

println("Mean innovation: $(@sprintf("%.4f", stats["mean_innovation"]))")
println("RMS innovation: $(@sprintf("%.4f", stats["rms_innovation"]))")

# Add convergence info
diag.iterations = 5
diag.converged = true
diag.cost_history = [100.0, 75.0, 55.0, 45.0, 42.0, 41.5]
diag.gradient_norms = [10.0, 7.5, 5.0, 3.0, 2.0, 1.5]

# QC stats
diag.qc_stats["ps"] = Dict(
    "total" => length(observations),
    "used" => stats["n_used"],
    "rejected" => stats["n_rejected"]
)

# Write diagnostic files
println("\nWriting diagnostic files...")
DW.write_all_diagnostics(diag)

# List generated files
println("\nGenerated files:")
for file in sort(readdir(output_dir))
    filepath = joinpath(output_dir, file)
    @printf("  %-40s %8d bytes\n", file, filesize(filepath))
end

println("\n✓ Test complete!")
println("="^80)
